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Abstract 

We describe a real-space approach to the calculation of the properties of an 
insulating crystal in an applied electric field, based on the iterative determi- 
nation of the Wannier functions (WF's) of the occupied bands. It has been 
recently shown that a knowledge of the occupied WF's allows the calcula- 
tion of the spontaneous (zero-field) electronic polarization. Building on these 
ideas, we describe a method for calculating the electronic polarization and di- 
electric constants of a material in non-zero field. The method is demonstrated 

for a one-dimensional tight-binding Hamiltonian. 
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Modern electronic band-structure and total-energy methods, such as those based on the 
local-density approximation (LDA) to density-functional theory are not well suited for 
studying the properties of materials in an applied electric field. These approaches typically 
rely on solving for the eigenfunctions of the effective one-electron Hamiltonian. Unfor- 
tunately, the electric field acts as a singular perturbation, so that the eigenstates of the 
Hamiltonian are no longer Bloch states, and the electron band structure is destroyed, even 
for arbitrarily small applied fields. Linear-response methods |@] are capable of computing 
derivatives of various quantities with respect to applied field, but cannot be used to study 
the electronic structure of a crystal in a non-zero electric field directly. Extensions of these 
methods beyond linear response [Q] consist of rather involved expressions which must be 
carefully handled in order to avoid divergences in the static limit. Formal studies of the 
structure of the Wannier-Stark ladder 0] are instructive, but are not of much practical help 
from a computational point of view. 

Recently, several groups have introduced new real-space approaches to the solution of the 
electronic structure problem. Based either on the locality of the real-space density matrix 
I^J^ or the use of a localized, Wannier-like representation of the occupied subspace 
these methods were motivated largely by the search for so-called "order-A^ methods" (for 
which the computational effort scales only linearly with system size A^). However, these 
methods also hold promise for application to the electric-field problem. The methods based 
on a Wannier-like representation |^,^ appear particularly promising in view of recent work 
showing that the electric polarization of a solid can be directly related to the centers of 
charge of the Wannier functions (WF) ||^. 

In this paper, we show that the real-space method of Mauri, Galli, and Car (MGC) 
can be developed naturally into a practical method for calculating the response of an 
insulator to an electric field. Both the spontaneous and induced electric polarization P are 
easily calculated, as are the perturbed charge density and polarization energy, and dielectric 
constants can be obtained by finite differences. The WF's are expanded in a local basis and 
truncated to zero beyond a real-space cutoff radius -Rc, these being the only approximations 
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involved. We apply the method to a one-dimensional (ID) tight-binding (TB) Hamiltonian 
and find that both the spontaneous polarization and dielectric constant converge quickly, 
with respect to i?c, to the values obtained by standard fc-space techniques While the 

method is demonstrated in a simple tight-binding context, we see no obvious obstacles to 
its implementation in a fully self-consistent ab-initio LDA calculation. 

We first summarize the MGC scheme using a density-matrix description. The exact 
density matrix of an insulating crystal is 

Pexact = I ^nk)(V'nk | , (l) 

nk 

where the {ipnk} are occupied Bloch eigenstates of the one-electron (e.g., LDA) Hamilto- 
nian. For insulators, we can define a unitary transformation from the {ipnk} to a set of 
localized (Wannier-like) functions {ipj} such that the occupied subspace is invariant under 
the transformation. In terms of these localized orbitals the density matrix is 

M 

Pexact = I] I I , (2) 

j 

where j runs over a set of localized orbitals that span the M-dimensional occupied subspace 
of Bloch eigenstates. The density matrix is spatially localized, i.e., p(r, r') as |r — r'| — >■ 
oo, the decay being exponential for an appropriate choice of phases of the Bloch functions 
0- 

In the MGC scheme, Eq. (^ is taken as an ansatz for a trial density matrix, p. Each ipj 
is free to vary within a real-space localization region (LR) of radius Rc, and is set to zero 
outside the LR. No orthonormalization constraint is imposed. The physical density matrix 
is then defined by 

P = pE(/-p)', (3) 

k=0 

where K is odd. The expectation value of any operator A is then given by tr[py4]. In the 
limit where p is idempotent we have p = p, as can be seen from Eq. (^. In general, the 
variarional p is not exactly idempotent, due to the approximations involved. The parameter 
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K controls the accuracy of the idempotency requirement [0]. Note that the set {V'j} in 
Eq. spans an M-dimensional subspace and that p spans the same subspace, since the 
function /x(a^) = 1 ~ (1 — x)^^^, corresponding to Eq. is such that /x(0)=0. Moreover, 
fxix) has its absolute maximum at x=l. From these properties it follows that: (i) when 
P = Pexact+0(5) then p = Pexact+C'((5^) << 1); (ii) the total energy, tr[pif], is variational 
provided that all eigenvalues of the Hamiltonian H are negative. 

In what follows, we restrict ourselves to K=l, so that p = 2p — p^. The variational 
total energy is given by £'[{■?/;}] = tr[p(2 — p)H], where H has been shifted to make all 
eigenvalues negative This functional is minimized with respect to the coefficients of 

the wave functions {ipj} in the given basis. In the limit Rc ^ oo, minimization of this 
functional yields a set of orthonormal orbitals that exactly span the occupied subspace of 
H, and consequently £'[{'?/'}]min = Eqs (the exact ground-state energy). For a finite Rc 
the resulting {ipj} are not exactly orthonormal and £'[{'?/'}]min > Eqs, as a result of the 
variational nature of the functional. 

We next propose a generalized formulation of the MGC scheme which is more natural 
and efficient in the case of periodic systems with small unit cells. In its original formulation, 
the MGC scheme is well suited to solving problems involving large supercells, where a 
minimization of a k-dependent energy functional, E''*^ [{■?/'}], is typically performed only at 
k=0. When applied to a periodic system with a small unit cell, one would have to minimize 
£'''[{'?/'}] independently on a mesh of k-points. In either case, the orbitals {ipj} are not truly 
localized, but rather are Bloch functions of the cell or supercell. Instead, we propose to let 
the wavefunctions be truly localized in the manner of Wannier functions (WF's), ipj = wi^n 
{I and n are cell and occupied-band indices respectively), and work directly with the original 
Hamiltonian. In this formulation, k plays no role. The electronic degrees of freedom are the 
coefficients of the WF's {ifo,n} in one cell. The periodicity of the system is now taken into 
account by introducing the periodic images of these WF's in neighboring cells, 

\Wi^n) = Ti |Wo,n) , (4) 
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where Ti is the translation operator corresponding to lattice vector R^. In terms of the 
density matrix, p = J2i,n \'^i,n){u!i^n\, the periodicity of the electronic state is expressed as 

p(r,r') = p(r + Rz,r' + Rz) . (5) 

The total energy per unit cell can be written explicitly as 

M M 

Eq [{w}] = 2 ^(wo,n I Ho I Wo,n) " H(^0,n | Wl^rn){wi,m \ Hq \ Wo,n) , (6) 

n=l n,m=l I 

where M is the number of occupied bands. Eq and Hq denote the total energy and the 
Hamiltonian at zero electric field. Since we are dealing with localized orbitals, only a finite 
number of orbitals will contribute to the second sum in Eq. (^. Minimization of Eq yields 
M approximate WF's, for a given choice of Re- Note that Eq. reduces to the original 
MGC scheme if the sum in the last term is restricted to /=0 only. 

It is thus seen that, when applied to a periodic system, the generalized MGC scheme 
is equivalent to the direct determination of its WF's. This brings us to the work of Refs. 
0, where it is shown that the electronic contribution to the spontaneous polarization of a 
periodic system can be written as Pg = — (e/fi) J2n'''^n, where r„ is the center of charge of 
the WF of band n in zero field, e is the electronic charge, and Q is the unit cell volume. In 
the present formalism, this becomes 

M M 



e 



n=l n,m=l I 



(7) 



At this point we have already completed the necessary steps for a real-space computation 
of the spontaneous polarization of a crystalline insulator, by combining Eq. (|^) with our 
modified MGC scheme. 

Now, we extend the relation between polarization and the centers of charge of the WF's 
to a periodic system in an electric field F. The Hamiltonian becomes 

H = HQ + eF-r, (8) 

Replacing Hq by H in Eq. @ leads to the total-energy functional 



E 



{w^} =Eo {w^} -eF-P, {w^} . (9) 



We retain Eq. and hence Eq. (^), and minimize this functional subject only to the 
constraint of localization of the field-dependent WF's {w^} to the LR, as before. Because 
of the locality of the WF's, the additional terms which enter Eq. through Eq. (|^) do not 
add appreciably to the computational effort. 

We emphasize that our solution does not correspond to the true ground state of the 
system. (There is no true ground state in finite field, as the energy can always be lowered by 
transferring charge from "valence-band" states in one region to lower-energy "conduction- 
band" states in a distant region.) We can think of our solution heuristically as the one 
which is generated from the zero-field state by adiabatically turning on F, and keeping the 
periodicity of the electronic state expressed in Eq. (|) . This should be very closely related to 
what is done when the field-dependent response of the crystal is measured experimentally. 
In order to measure "static" properties, the field must be turned on (or allowed to oscillate) 
on a time scale that is slow compared to usual electronic processes, but fast compared to 
the characteristic electronic tunneling rate at the maximum field encountered. Thus, the 
experimental object of study is also really an excited state (more properly, a very narrow 
resonance) of the finite-field Hamiltonian. Our solution has a very similar interpretation. 

One might object that the existence of a functional -Eoll""^}]; as in Eq- (H), relies on 
the identification of bands. In the presence of an electric field, where all gaps disappear 

the existence of bands is not obvious. Nevertheless, Nenciu [|] has shown that one can 
define a sequence of periodic Hamiltonians {Hq}, constructed from H by projecting out the 
non-periodic part of F ■ r, such that the subspace of occupied states of a given Hg reduces to 
that of Hq in the limit F — 0. It is argued that the "bands" defined by the Hamiltonian Hq 
provide an increasingly accurate description of the finite-field electronic state as the integer 
index q is increased. (Eventually, as q gets too large, the behavior diverges, in the usual 
manner of asymptotic perturbation theory. In other words, the radius of convergence Fq 
tends to zero as g — * oo.) At least at small F, our {w^} are presumably very similar to the 
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Wannier functions that would be constructed from the bands of Hg, and thus should give a 
good description of the experimental electronic state of the system. We return to this point 
below. 

Before turning to our tests and results, we mention one technicality. The minimization 
of the functional of Eq. can be performed directly at fixed F using steepest-descent 
or conjugate- gradient techniques. This appears to work quite well at weak fields, but can 
become unstable for strong fields. Alternatively, one can perform the minimization with a 
constraint of fixed (r) (i.e., fixed P), treating F as an adjustable Lagrange multiplier. Since 
(H) — F - (r) defines a Lagrange transformation from (H) to (Hq), it follows that V(r) (-f^o) = 
— F. In this way, the function F(P) can be mapped out, and inverted numerically to give 
P(F). The latter approach is more appropriate for investigating the strong- field behavior of 
the solutions. 

We apply our scheme to a ID tight-binding 3-band Hamiltonian in which each unit cell 
consists of three atoms with one orbital per atom, 

H (a) = ^ jcj (a) c^jCj + t c]cj+i + h.c. | , (10) 

j 

with the site energy given by e^m+kici) = Acos{a — Pk)- Here m is the cell index, k = 
{ — 1, 0, 1} is the site index, and f3k = 27rA;/3. This is a simple model of a sliding commensurate 
charge- density wave which slides by one period as the parameter a evolves through 27r. For 
our tests we set e = 1 and t = A = —1, and use x = with xj = j/3 for the 

TB position operator. We discuss the results obtained with only the lowest band filled; the 
discussion applies equally to the case of two filled bands (the 3-filled-band case is trivial). 

We first consider the TB Hamiltonian with no external field, and calculate the sponta- 
neous polarization as a function of a, comparing with the results obtained by the method 
of Refs. which we take as exact for the sake of comparison. The results are shown in 
Fig. |l](a). In the inset we show the exact polarization as a function of a in the interval 
[0°, 120°]; also plotted are the results for |Pcxact(«) — -P(«)| in the same interval, for Rc = 1.5 
(9 sites within LR) and Rc = 2.5 (15 sites). The LR was kept centered at the origin for all 



values of a. Convergence is already very good for i?c = 1-5 with a maximum percentage 
error of ~1.5%, dropping to ~0.5% for R^. = 2.5. 

We now apply an external electric field to the system. We minimize Eq. (|^) for six fixed 
values of F between 0.01 and 0.06. In this field region the polarization P is linear with 
F, to a very good approximation. In Fig. |I|(b) we show the linear dielectric susceptibility 
X as a function of a for Rc = 1.5, Rc = 2.5 and Rc = 3.5 (21 sites). Also shown are the 
linear-response results we obtained using the method of Ref. 0. x converges less rapidly 
than P, but the maximum percentage error for a G [0°, 120°] is already ~3.0% for Rc = 3.5. 
Convergence is systematically worse around a = 60° where the WF's are least localized, due 
to the fact that the gap between the two lowest bands reaches its minimum value of 0.814 
at a = 60°. 

Next, we minimize Eq. keeping (x) fixed. For a given value of Rc we explore the 
interval [—Rc, Rc] of possible values of (x). Fig. shows Eq = {Hq) as a function of P = (x) 
for a = and Rc = 1.5, 2.5, and 3.5. For unconstrained WF's we expect Eq to be periodic 
with P. Our results reproduce that behavior remarkably well, except near the boundaries 
of the LR where the variational solution becomes poor. Note also that for a given value 
of P, the value of Eq diminishes with increasing Rc, reflecting the increasing quality of the 
variational solution as Rc increases. 

We point out that it should be straightforward to obtain the higher-order dielectric con- 
stants from the Eq(P) curve, by performing careful finite difference calculations. Moreover, 
the generalization of our approach to an ab initio LDA calculation should be easily imple- 
mented. A localized-orbital basis would be ideally suited, although a plane-wave basis might 
also be used ^j^. As in the original MGC scheme, self-consistency can be included in a 
straightforward manner, since the density n{r) = p(r, r) remains periodic. Thus the Hartree 
and exchange-correlations terms can be computed as usual and do not contribute to the 
non-periodic part of the self-consistent potential. 

A final word of caution is in order. It should not be imagined that there is a well-defined 
curve Eo{P), periodic in P, which can be obtained by taking the limit -Rc — ^ oo of our 
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procedure. On the contrary, as the LR grows very large, it becomes possible to construct 
a solution for w having arbitrary [x) (i.e., arbitrary P), and energy arbitrarily close to 
Eq{F = 0). This can be done by starting with the zero-field WF of the occupied band and 
admixing a small amplitude (of order /~^/^) of a zero-field unoccupied-band WF at a distance 
/; its energy approaches Eq{F = 0) as / — oo. Thus, we have the pathological situation that 
Eq{P) becomes perfectly fiat in the limit Rc —* oo. (When working at fixed F, this pathology 
shows up in the form of a growing number of false local minima of the functional of Eq. (^ 
as -Rc oo.) However, the Taylor coefficients of Eq{P) (expanded about the minimum) are 
well-behaved in the limit Rc —>■ oo, even as its radius of convergence is decreasing to zero. 
Underlying this behavior is the asymptotic nature of the expansion, which is also the case for 
the Nenciu construction of the "polarized" sub-spaces. In both cases, full convergence as 
Rc oo or q ^ oo is obtained only in the limit F — 0. We do not claim that our proposed 
method has superior convergence or gives more physical results than that of Ref. |^ as F 
gets large. But both approaches must have the same small-field behavior, and hence yield 
the correct perturbation coefficients (e.g., linear and nonlinear dielectric constants). The 
proposed method also has the advantages of being computationally tractable and convenient 
to implement. 

In conclusion, we propose a method for calculating the response of an insulator to an 
applied electric field based on a Wannier-function-like representation of the electronic or- 
bitals. In this approach the spontaneous polarization, the perturbed charge density and the 
polarization energy are easily obtained, and dielectric constants can be calculated by finite 
differences. The method is variational, and therefore is well suited to solution by iterative 
techniques such as conjugate gradients. The computational effort scales only linearly with 
system size and the method becomes exact as the cutoff radius used to truncate the Wannier 
functions is increased. The method is demonstrated in a simple tight-binding context, but 
is also well suited to implementation in a fully self-consistent ab-initio LDA calculation. 
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FIGURES 

FIG. 1. Polarization P{a) and linear dielectric susceptibility x(a) for the tight-binding model 
of the text, (a) \P{a) — -Pexact(")| for Rc = 2.5 (dashed) and Rc = 1.5 (dotted line). Inset: 
-PexactCo;)- (b) x(q!) from linear response (solid line) and from current method with Rc = 3.5 (long 
dashed), Rc = 2.5 (dashed), and Rc = 1.5 (dotted line). 

FIG. 2. (Hq) as a function of (x), for a = 0. Solid line, Rc = 3.5; dashed line, Rc = 2.5; 
dotted line, Rc = 1.5. 
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